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We apply phase space analysis to inhomogeneous cosmological model given by Lemaitre-Tolman 
model. We describe some general conditions required to interpret the model stable enough and, in 
the present paper, apply them to two special cases: dust filled homogeneous model with and without 
cosmological constant. We find that such stability explaining all present astrophysical observations 
can not be achieved due to instabilities in phase space. This hints that non-homogeneous models are 
not likely to be physically viable, although any conclusive analysis requires more realistic modeling 
of non-homogeneous universe. 

I. INTRODUCTION 

Supernovae observations^! made just before the break of the millennium implies that the universe appears to be 
expanding at an increasing rate. A bit later made observations of cosmic microwave background (CMB) radiation 3 
supports the conclusions made out of the supernovae observations. The most popular ways of explaining these 
observations are with models based on Friedmann-Lemaitre-Robertson- Walker (FLRW) metric, which is based on 
general relativity and the principles of cosmological and Copernican. The cosmological principle merely states that 
the universe is spatially homogeneous everywhere, whereas the Copernican principle states that there are no preferred 
points in the universe. FLRW metric and models based on it are extensively presented in the literature of cosmology 
(see for example^). 

Even though the FLRW based models fit well inside the frame pr ovid ed by our observations, it is also for l ong 
known to suffer some problems, for example the fine tuning problerrl^ and the cosmological constant problerrFEJ 
One of the strengths of the FLRW based models is simplicity due to varies approximations, although, that can also be 
counted as a weakness. It is also questionable if all approximations are made acceptably, as is pointed out by Shirokov 
and Fished, where is questioned if the homogeneity approximation should be done to the Einstein tensor G^, rather 
than to the metric g^ v , since in general < Gn V {g^) >^ G^{< g^ >). However, the problem is more complex than 
this. As pointed out by Shirokov and Fished, Einstein equations are no longer tensor equations after averaging in the 
sense, that they can not be changed e.g. from covariant form to contravariant form with metric tensor without altering 
the equations. In this perspective it seems, that only tensors rank and scalars have well validated averages. For this 
kind of approach see BuchertP, where he transforms the Einstein equations into scalar equations before averaging. It 
is very much possible that the problems FLRW based models suffer are due to approximations, and solely the work 
of Shirokov, Fisher and Buchert implies that the first approximation to investigate more is the homogeneity. Such 
a model was first introduced by Lemaitre^, and later on studied by TolmarP^ anc l its called Lemaitre-Tolman (LT) 
model. For further developments of the model se c! 11 * 12 *. The LT model have not yet been studied as widely as FLRW 
based models, hence all the problems it suffers have probably not yet been discovered, but it have already shown its 
power by overcoming some of the problems of the FLRW based models. For example, Mattsson has shown that the 
LT model can explain the main cosmological observations without dark energyfSl 

Actually, FLRW based models are often presented including early times inflation, which solves the fine tuning 
problem. This solution can not be generalized to LT model a prior, because the evolution of the universe in LT 
model is dependent in coordinate distance, which would make the inflation occur differently in separate locations, and 
the consequences of this are unknown. However, inflation is not the only possible explanation for the homogeneous 
tendency of cosmological observations. We explore the possibility that the structure of the equations governing the 
evolution of the universe is such, that it has a inbuilt property to make everything appear as observed. We study the 
existence of this property by using phase space analysis. The aim is to find restrictions to viable and stable solutions 
of the differential equations governing the universe in the LT model. Viable solutions mean here that the solutions 
are consistent with the observations. Stable solutions include such solutions, which are attracted towards the universe 
we observe; therefore these kind of solutions do not need fine tuning of initial state. Especially interesting cases are 
where the dark energy is absent. 

It is also interesting to see, if our results offers insight to the homogeneity approximation. By that we mean, if all 
the viable and stable LT models are approximately homogeneous. This subject however is not going to be important 
in this paper, but it is merely pointed out as a possibility what more can our results offer. 

The focus of the present paper is to introduce a novel method to use stability analysis and its general features with 
the Lemaitre-Tolman model, which is done in section II. In section IIA, homogeneous cases in general are applied to 
the methods found out. In sections IIB and IIC pressure free, flat and homogeneous universe is investigated, in the 
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cases of dust filled and dust and dark energy filled universes, using the phase space analysis. Finally in section III the 
results are discussed. Realistic applications where viable and stable inhomogeneous models are to be determined are 
left to forthcoming publications. 



II. LEMAITRE-TOLMAN MODEL AND PHASE SPACE ANALYSIS 

The LT modeP^ describes a dust filled inhomegeneous but isotropic universe which energy momentum tensor reads 
as = pu^u^W^Here p — p(t, r) is matter density and is the local four- velocity. As the coordinates are assumed 
to be comoving, the four- velocity is simply = 8^. The standard synchronous gauge metric in the LT model is given 

ds 2 = _ dt 2 + R r d ^ + R{t] r) 2 l d6 2 + sin 2 q d ^ (1) 

l + /(r) 

where the subscript r denotes derivative with respect to radial coordinate, R r = dR(t,r)/dr, and /(r) > — 1. In this 
prescription the evolution of universe is built in to the local scale factor R whereas function / controls the overall 
radial dilatation. Including the cosmological constant A into the Einstein equations, after some integrations, one 
obtains the relevant differential equations as 

R 2 = 2 -f + f + ±R 2 , (2) 

and 

KP= R-W> (3) 

where M = M(r) is a arbitrary function, k = 87rG/c 4 , and G is the Newton's constant of gravity. 

In general, all the quantities (or functions) in Eqs. ([2| and ^ are unknown. However, if they are presented as 
quantities dependent only on redshift, they can be received from cosmological observations, or be derived from the 
obsevations, e.g. if R(z) is known we get R z (z) as its derivative. In the LT models the redshift equation reads as 

efe _ (1 + z)Rrt(r, t) 

dr VTTJ ' ( ' 

and the angular diameter distance d a is d a (z) = R(t(z),r(z))^HIB 

The relation between t and r can be given as follows. The path of radial light ray is given by the radial null geodesic, 
where ds 2 = d6 2 = dijj 2 = 0. Using metric |l]) it is 

dt = ±^=dr, (5) 

where the signs correspond an incoming (— ) and an outgoing (+) light rays. It is possible^ to choose radial coordinate 
r, i.e. use remaining gauge freedom so, that for incoming ray dr = —dt. So, along the ray 

t = t{r) = t - r, (6) 

where t = to refers to present time. For any quantity, a hat " denotes that it is evaluated along (incoming) light ray, 
e.g. R = R(t(r),r). Thus, the gauge condition for Eq. ^ for incoming light can be written as 

1 = . (7) 

On the radial null geodesic it can be showri^ that the relation between the matter density p(t, r) and the number 
density of light sources in redshift distance 23 n = n(z) is given by 

- m dz 

P =h 2 ^ (8) 
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whereu = n(z) is the mean mass per source at given redshift distance. After some manipulation, using Eqs. j3j, 
ffify, Q, and ([8]) the redshift equation Q can be cast in the forrrP 

J(z)z r + K(z)z 2 +L(z)z rr = 0, (9) 

where we have defined^ 

J(z) = k(1 + z)/in, 

K(z) = 2R[R z + (l + z)R zz ], (10) 
L(z) = 2(1 + z)RR z . 

All the quantities in Equations J, K, and L are only dependent on redshift z, and they are all measurable or they can 
be derived from measurable quantities. Moreover, because it is evident that the quantities dependent only on redshift 
are evaluated along light ray, every function or quantity dependent only on z is presented without a hat. Especially 
now R = R{z) = R(t(r(z)),r(z)), and hence R z = dR(z)/dz = dR(z)/dz and R zz = d 2 R{z)/dz 2 = d 2 R(z)/dz 2 . Eq. 
([9]) is a second order non-linear differential equation, from which can be solved z as a function of r. Howeve r, tha t is 
not the interest here, but rather investigating the stability of different solutions. We use phase space analysiiP^^D f or 
investigation. Here is chosen z r = y to give 

_ J(z)+K(z)y 

~ y W) (11) 
= y- 

Note, that we have given the redshift in terms determined comoving coordinate r, which is related along the ray to 
time as r — to — t. This means that we can use practically interchangeably the two parameter r or t. 

Analyzing cosmological data will give us functions J, and L explicitly with respect to z, but even now when the 
explicit forms are unknown, general remarks can be done on what kind of systems Eq. ([£]) can describe. The phase 
plane can be divided into sections each having its own characteristics. Some of the properties can be read out from 
general properties of the equations and represented as a flow plot in zy-plane. 

The y-axis represents location of observations (Earth) at present time; the redshift is there zero. On the right hand 
side of the y-axis is past or distant objects and the left hand side of the y-axis can be interpreted as future. On 
the z-axis z r = z t — 0, hence it represents apparently static universe. Above the z-axis z r > 0, the area represents 
apparently (and locally) expanding space, and below the z-axis z r < the apparently contracting one. Therefore, 
due to the observations, we are mostly interested of the first quarter of a phase space plane, past of an expanding 
universe; however, we do not want to exclude other parts of the phase plane per se, but we give the first quarter of 
the phase space plane most of our interest. 



Looking to the latter of the Eqs. ( 11 ) one sees, that whenever y > the flow arrows points right, and when y < 



they points left. The curves where y r or z r is zero are called nullclines and the points where y r and z r is zero are 



called fixed points. In the system ( 11 ) all the points on curve y = are fixed points, which means that the curve y = 



is fixed. The curves where L(z) = needs special attention as the system ( 11 ) is not well defined there. The physical 
perspective also gives more restrictions. The form of function J reveals that it is positive implying that the curve 
J{z) + K(z)y = can not cross the z-axis. The only exceptions are at the origin of the spherical symmetry (r = 0), 
where it is required to have R(t, 0) = for all t to avoid point mass and curvature singularity at r — 0, and possibly 
at the Big Bang or the Big Crunch!^ Requirement R(t, 0) = in gauge (|6| is R(to, 0) = 0, which is compatible with 
R{z) being zero at z — as R(z) is the angular diameter distance. Because J is always positive (neglecting the special 
points discussed above) and K and L seem to able to change their signs, combining these functions with different 



allowed signs there is essentially four different types of situations (in the case of system ( 11 )) that can appear on the 
phase portraits, presented in Fiq. [lj The nature of the function R(z) implies that L is zero only if R z (z) is zero 
(neglecting again the special points discussed above) and because R z {z) — R r {r)r z , is then either R r {r) or r z zero. 
The latter case can approach to zero when y — > 00, and R r {r) — occurs at the apparent horizon^. 



Consider now the situation where functions J, K, and L are formed from observational data. Eq. (11) is now 
given explicitly, so phase plane can be drawn and the best fit curve compared to observations can be fitted. It is 
assumable that the best fit curve coincides closely to the best fit to the observations in the isotropic and homogeneous 
models, which means that the best fit curve is not very "lumpy" and is approximately monotonic. Now, consider each 
observable separately. Each of them have their own trajectory, which are unknown to us, because we do not have 
enough data from each observable. But if we had enough data from each observable, we could draw their trajectories 
on the same phase plane that the best fit curve is on. So, each observable have their own trajectory, but still they 
seem to sit on a approximately monotonic curve, the best fit curve. In our scheme this appears to happen only, if 
the flow arrows are pointing towards the best fit curve. This would ensure that for large number of initial values, 
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FIG. 1: Four different situations depending on the signs of K and L that can occur on a phase plane with the system Functions J, 

K, and L are here chosen to be +1 or —1 for simplicity. The dashed curve is the nullcline J + Ky = 0. 



observables end up (after long enough time period) near by to the best fit curve. This issue, however, is not that 
simple and it will be discussed more in section Conclusions and discussion. 

At simplest the interest in phase plane analysis is concentrated into the existence and properties of the fixed points. 
In this case however, observations suggests that the best fit curve should be attractive. In the system (11 1 only z-axis 
can be referred as a attractive curve, but its nature as apparent static universe is not what is observed. However, 
nullclines can be thought of to be attractive like, since flow arrows can point towards them as in e.g. the low-right 
case in Fig. [I] In our system the only physically interesting possible attractive like nullcline is J(z) + K(z)y = 0. 
However, the exact identification between the nullcline and the best fit curve can not be done, because then both Eqs. 
J + Ky = and ^ should be satisfied simultaneously causing either the solution y{z) be linear or L{z) = 0, but it 
is not even necessary as long as the identification can satisfy observations inside their inaccuracy limits. 

The properties of the system (11) discussed above implies that the curve J + Ky = can be attractive in the 
first quarter of the phase plane only if J is positive and K and L are negative. Even though these boundaries are 
necessary, they are far from sufficient for the following reason. Consider a situations where < J, and K, L < and 
the value of the slope of the curve J + Ky = 0, —J/K, is positive. On the curve J + Ky = is always y r = 0, thus 
flow arrows on the curve are horizontal and pointing right, i.e., the slopes of the flow arrows y z — y r /z r are zero. 
Now, approaching the curve from below vertically makes y z to approach zero. At some distance to J + Ky = before 
y z reaches to zero, it becomes smaller than — J/K. This means that the solution is no longer approaching the curve 
J + Ky = 0, because flow arrows determine the slopes and the progressing directions of the solutions in each point. 
Similar situation occur, if at fixed z we approach the curve J + Ky = from above and — J/K < 0. This is also why 
nullclines are rather attractive like than attractive. As it follows, we need more sufficient methods to measure the 
attractivity of the nullcline J + Ky = 0. Especially, we need a method to measure the distance from the nullcline 
where it stops acting as an attractor. In the present paper for this is used the following method. 

Let us consider our system on area z, y > with positive J and negative K and L. The nullcline J(z) + K(z)y = 
is attractive like, if flow arrows around it points towards it. Now below the nullcline flow arrows point up-right and 
above the curve they point down- right. The slope of the nullcline is positive and all the flow arrows above it points 
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towards it, but the flow arrows below it points towards it only if the slopes of the flow arrows y z {z) at some given 
z are greater than the slope of the nullcline — J/K at the same z. Let S y > be the (vertical) distance from the 
nullcline J{z) + K(z)y = 0. The slope of the flow arrow at given z and distance 6 V below the nullcline is: 



(V - 5 v)r 



y=-J(z)/K(z) 

J(z)+K(z)(y+S y ) 
L(z) 



y-S y 
J(z)+K(z)(y-S y ) 



y=-J(z)/K(z) 



K{z)6 y 
L(z) 

K(z) 



L(z) 

■ J(z) + K(z)y 
L(z) 



y=-J(z)/K(z) 



y=-J(z)/K(z) 



(12) 



From Eq. (12) one can explicitly see, that the size of the vertical step S y taken from the nullcline effects to the 



attractiveness in a linear fashion: the larger the step is, more attracitvc like the nullcline seems. This means that S y 
takes a role of a parameter comparable to observational inaccuracies. Hence, the slopes of the flow arrows y z {z) are 
greater than the slope of the nullcline —J/K at at some given z, if 



L{z) v>Vz K(z) + K{zY z[h 
The above inequality is the restriction we use to study the attractiveness of the nullcline J + Ky = 0. 



(13) 



A. Isotropic and homogeneous models 



To give more concrete touch of our prescription, we check how the method works with isotropic and homogeneous 
models. Note, that FLRW metric can not be used here, because it is not compatible with our gauge choice dt — —dr. 
This can be explicitly seen by comparing the standard FLRW metric relation between r and z given bj^ 



_ 2zn + (2n - 4) (Vz^o + 1-1) 
T ~ a H (z + l)n 2 ' 

with a = do/(l + z ) t° the relation derived from Eq. (M), a = vT— kr 2 giving 



r = ±i 



k k(l + z) 2 ' 



(14) 



(15) 



Clearly Eqs. (14 1 and (151 are not equivalent, except in some special cases. However, it is necessary to write isotropic 



and homogenous space-time in Robertson- Walker coordinates, and we can proceed by other means. 
In spherically symmetric homogeneous space-time the metric can always be given as^ 

fc(u • du) 2 
1 - fcu 2 



ds z = g{v)dv 2 + f(v) du z + 



(16) 



where v and u = (ui, ^2,1*3) are the coordinates, g(v) is a negative and f(v) is a positive function of v, and k is 

spatial curvature and can be chosen to be 1, 0, or -1. For our purposes it is convenient to define new coordinates t, 
r, 9, and ip by 

yj-g{v)dv = dt, 

ui — O(r) sin 9 cost/?, 

U2 = 0(r) sin sin 95, 

u 3 = 6(r)cos0. (17) 
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Then we have 



ds 2 



-dt 2 + a 2 (t) 



ft 2 



e 



i-fce 2 



dr z + e 2 dn 2 



(18) 



where a(t) = J f(v) and dQ 2 = d9 2 + sin 2 9 dip 2 . With the above metric definition the field equations take the normal 
Fricdmann form and therefore can be written as 



a t =H a /^Of (a/ao)-^ 1 ^), 



(19) 



where is the present value of the energy density of different energy forms. For dust, dark energy, and spatial 
curvature, the Wi takes values 0, —1, and —1/3 respectively. Because the LT model does not take pressure into 
account, for our purposes it is unnecessary to include relativistic matter here either, even though generally it could 
be done. 



To be able to use the Eqs. (11 ), we need to specify quantities p, n and R. According to Eq. (18]), p and n can be 



given with R, p, and z r , and because in metric (18) R = a<d, we need to find out presentations for quantities a, 0, p, 
and z r with respect to z. 

In LT model, the energy momentum tensor T^ v includes only dust, hence 

p = p^(a /a)\ 

where p^ — fl ffi ^ fc- , and £1$ is the present dust (or non-rclativistic matter) density of the universe. 
With metric ( 18 i equations |7|) and Q reduces to be 



and 



~ Vi - fee 2 

dz (1 + z)a t @ r 
dr ~ Vl-fce 2 ' 

Combining these the usual Robertson- Walker relation 

1 a 



1 + z ao 

is reproduced and as r — to — t is z r = aoa t /a 2 , we can write 



(20) 



Zr = H {z + l) f£n\ 0) (z + l)3d+ 



(21) 



(22) 



(23) 



(24) 



In the given gauge, the function O can be calculated directly from Eq. (21 ), which can also be written as 



dQ 



Vi - fee 2 



1 + z 1 

o a o z r 



-dz, 



where we have used the relation 



r dr 



o a Jo a dz 



1 dr 



-dz 



1 + 2 1 

a z r 



-dz . 



(25) 



(26) 



0(r = 0) = is determined by the requirement R(z = 0) = 0. The solution for curvature cases can now be integrated 
out. We find 



O = sinh 



1 + 2 1 

a z r 



-dz 



-I. 



(27) 



°0 z r 



k = 0. 



(28) 
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= sin 



1 + z 1 



dz 



k = l. 



(29) 



In the following sections some numerical calculations are executed. To carry out the numerics we have chosen 
lOOkm/s/Mpc = 1, which makes all the physical quantities dimensionless. For subsequent calculations also it is 
convenient to define a function 



1 

Ho 



J z (z) , J(z) 



K{z) K{z)' 



K z {z) 



L{z) 

K( z y 



(30) 



which in homogeneous models is g — g(z;tt°,ao,Ho), i.e., dependent on variable z and parameters £1®, op, a nd Hq. 
Even further, it is easy to see, that in flat homogeneous cases g = g(z; fi?), thus the boundary condition (13) can be 
written as 



%->gM). 

-no 



(31) 



The fact that g(z; Q®) do not include parameters ap, Hq, or k, considerably simplifies the boundary condition. 





FIG. 2: Phase space portrait of the dust filled homogeneous universe (left hand side plane), when Ho = 0.70, Q m = b and Q\ = 0, and 
of the dust and dark energy filled homogeneous universe (right hand side plane), when Ho = 0.700, Q m = 0.279, and Q\ = 0.721. In both 
planes: thin purple curve is J + K y = 0, thin black curve is L = 0, and red, green and blue thick curves are solutions of the system 
with Hq values 0.69, 0.70, and 0.71 respectively. 



B. Dust filled universe 



Let us study stability of the flat dust-filled homogeneous universe, which is often used approximation to investigate 
the evolution of the universe during matter dominated era. From Eq. ( 24 ) we obtain redshift relation 



where f2 



(0) 



z r = H Q (z + 1)^^(1 + 2)3, 
1 is set from now on. The pair of differential equations dill) is now 



V 


[y 


(9z - V(z + l) 3 + 9^ 


- 


- 6H (z + l) 2 (z{z + 2) - y/(z + 1)3 + l) ) 


2(z + l) 


[iz - 2^/(z + l) 3 + 3) 





(32) 



(33) 



which is defined for non-negative z except when L(z) = at z = 5/4, which corresponds to R z — as can be seen 
from (10). Curve J + Ky — of the system (331 can be attractive like on the first quarter of the phase plane when 



L and K are negative and J is positive, which occurs at range 5/4 < z < 65/16, assuming H > 0. A situation of 
this kind is illustrated in Fig. [2j The case in the figure does not show any signs of attractivity: it seems that at the 



range 5/4 < z < 65/16 none of the flow arrows below the curve J + Ly = would lead any solutions towards it. This 
deduction is strengthen by (31), which (by assuming Hq > 0) in this case reduces to be: 




l) 5 (7339 - 8z (64z 2 - 516z + 963)) 2 8(z(112z - 883) + 1192) (s + 1) 



(65 - 16z) e 



(16z - 65) 3 



< 



"y 
H : 



(34) 



where the subscript m marks that this is the explicit form of the function g(z; ftj) in the dust filled case. At the range 
5/4 < z < 65/16 the function g m (z) is monotonically increasing from zero to infinity, which suggest that the nullcline 
,/ + Ky = of the system (33) can hardly be attractive like throughout the gap 5/4 < z < 65/16. 

The left hand side portrait in Fig. [2] illustrates what happens close by z = 5/4, where the system (33) is not defined. 



Even a slightest variation from today's observed value of Hq seems to lead to a very different kind of universe. It is 
notable, that numerical calculations do not give the same answers for Eqs. (33) and (32) even with Hq = 0.7, since 
the solution of Eqs. (33) can not cross the point z = 5/4. The area where J is positive, and L and K are negative, 
does not seem attractive at all, and only very close to the singularity it might be attractive like according to g m (z). 



A:, 
6000 h 
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FIG. 3: Both functions, g\ m and z°*, can be interpreted as a difference between two z r values, and are presented as such in a (Az r ,z)- 
plane. In both planes the solid curve is g\ m and the dashed curve is z° % . On the left hand side are ja„ and z°' drawn throughout the 
values of z where the nullcline can be attractive like, and on the right hand side are the same curves drawn at 1.63 < z < 2.00. On the 
left hand side plot the scale of Az r is such, that in practice the dashed curve is indistinguishable from the z-axis. From the right hand 
side plot one can see that at what range of z values is Sy/Hg smaller than observational inaccuracies. 



C. Dark energy and dust filled universe 



Next we assume universe to be homogeneous, flat and consists of dust and dark energy, thus from Eq. (24) we 
obtain 

Z r = H (Z + l)y/(z + l)3ft m + fi A . (35) 

Quantity R(z) takes now a very inelegant form, which also is the case in every function including R(z), hence none 
of the functions including R(z) is explicitly presented in this subsection. Also, R(z) includes an elliptic integral of 
the first kind, which makes accurate algebraic analysis difficult. However, numerical methods are sufficient enough in 
this case. 

Observations restrict the parameter values of this model: six-parameter ACDM model fit to WMAP nine-year data 
gives ft A = 0.721 ±0.025, ft m = 0.279 ±0.025, and H = 0.700 ±0.022^1. From Eq. ((35} can be seen, that the effect of 



dark energy with ft m = 0.279, ftA = 0.721 is less than one percent when z > 5.37, which is less than the inaccuracies 
of observed ft m , ftA, and Hq values. Hence, at 5.37 < z we approximate this model model to have ftA = 0. This 
approximation and numerical analysis at < z < 5.37 reveals, that L{z) — only at z 1.63, and L and K are 



negative and J is positive approximately at 1.63 < z < 5.18. In this case the boundary condition (31) reduces to be 
of the form 

77- > »Am(z; ftA, ftm), (36) 
-no 

where the subscript Am marks that this is the function g(z; ftj) in the dust and dark energy filled case. The explicit 
form of the function g\ m (z; ftA, ftm) is neither elegant nor necessary to show here, thus it is not presented here. At 
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1.63 < z < 5.18 the function g\ m {z) is monotonically increasing, and it is increasing in a very rapid fashion (as can 
be seen from Fig. p$b. IrPSl are the values of and £l m given with error margin ±0.025, so function 

m \z r (z; fi A = 0.696; fi ro = 0.304) -z r (z;Q A = 0.746; O m = 0.254) | 
z,, := (37) 

represents the observational inaccuracy independent of parameter Hq. From the portrait on the right hand side of Fig. 
[3] can be seen, that g\ m (z) < z° l at 1.63 < z < 1.70, which therefore is the gap where this model can be attractive 
like. 

The time period corresponding z values 1.63 — 1.70 takes place about 9.70 — 9.85 Gyr ago. Even though about 
0.15 Gyrs could be enough time for universe to evolve into similar state everywhere, thus possibly explaining the 
observations done at some range of z values, it certainly can not explain the homogeneity of the observations at 
larger redshift values, e.g. over 10 Gyr old galaxies or CMB. The system with values f2 m = 0.279, Ha = 0.721, and 
Hq = 0.700 is illustrated on the right hand side on Fig. [2j The phase portrait looks similar to that of the dust filled 
case. 



III. CONCLUSIONS AND DISCUSSION 



Phase space analysis is used to find out if there is isotropic dust source models describing the universe without fine 
tuning the initial state. In general with inhomogeneous models, it is not always reasonable to use term initial state 
to refer only to the state of the system at Big Bang time, but sometimes also later times. This is due to observations: 
there is no observations of each location at all time during the entire history of the universe. The Big Bang might of 
occurred at different times in different regions, and the regions might of evolved completely differently compared to 
each other. Hence, here the term is not fixed to some particular time or place, but is rather determined for each case 
separately, whatever is most convenient for the situation. Often we use it to refer to the time relatively near in the 
past of each observable, but there are some exceptions, such as homogeneous case, which are discussed more later. 

When cosmological observations are plotted on the (z r , z)-plane, the curve that is the best fit to those plots is 
the curve we refer as the best fit curve. An cosmological observable to be observed in the vicinity of the best fit 
curve today, does not demand that the observable have to evolve towards the curve throughout its history; it only 
needs to be there when observed. A single observable could in principle move repeatedly towards and away from the 
curve during its entire history, as long as it locates close enough the curve when we observe it. Though the amount 
of observables basically eliminates this kind of excessive behavior, but a reasonable situation would be where each 
observable at some point approach the best fit curve and stay close enough before observed. 

The interpretation of the dynamical equations simplifies in the case of homogeneous models, as the initial state 
is same everywhere and the Big Bang time occurred simultaneous everywhere; none of the quantities in equations 
governing the evolution are dependent of r, and the evolution have been similar always and everywhere. Thus, the 
minimum requirement for a FLRW metric based on model to be viable and stable enough is, that there has been a 
period where it have evolved towards the best fit curve and located inside boundaries given by observation inaccuracies 
since. This is the minimum requirement in the sense, that there was at least one early period where the best fit curve 
was attractive. The given example of the dust and dark matter filled case have an era where redshift is approaching 
the best fit curve: about 9.70-9.85 Gyrs ago. Hence the model could explain almost the last 10 GyrsPS Thus the 
model can not explain oldest observations. This is not surprising: as it is well known, without inflation homogeneous 
models have a fine tuning problem of matter and energy densityP^ In the light of this paper, the stability is not a 
inbuilt feature at least in the cases studied in this paper, and can not replace the widely excepted solution of early 
time inflation for the fine tuning problem. 

Analogous results are expected for almost homogeneous cases, where initial state does not vary too much with 
respect to r and the evolution can be though of almost similar everywhere: observations can be explained without fine 
tuning the initial values, if at some point in the history before CMB there was a period where the universe evolved 
towards the best fit curve (and have been close enough it since). Thus we can conclude, that if the observed universe 
have not had a period in its history before CMB where the boundary condition we gave is satisfied, LT model can 
not explain its homogeneous nature without fine tuning. 

The redshift values where the system is not defined, singular limits, causes difficulties for analysis. In this paper 
singular limits are not studied, and it is left to forthcoming publications. Nevertheless, one note of singular limits 
should be brought out here. In fig. [2] is illustrated what happens close by z = 5/4, where the system is not defined. 
Even a slightest variation from today's observed value of Hq seams to lead a solution into very different direction at 
the singular limit. On the other hand, if we would of drawn the phase plane with e.g. Hq = 0.71, the solution with 
H = 0.71 would of approached towards point ~ (1.25, 5) and the solutions with Hq = 0.70 and H = 0.69 would of 
approached towards the z-axis. This may be a sign, for example, of an instability of the model or of a break down 
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of the used method at this z value and its close neighborhood. As singular limits in general, also this phenomena is 
planned to be investigated thoroughly later. 

The overall view our analysis cast on inhomogeneous models is not very promising. The best fit curve that should 
be attractive like is approximately given with Eq. (35), thus increasing at rate z 5 ^ 2 with large z values, and rapid 
increase of the curve makes it more unlikely to be attractive like. The situation seems to get worse if pressure is 
taken into account, since the increasing rate is then proportional to z 3 with large redshift. However, it is possible for 
observationally acceptable inhomogeneous (and therefore also homogeneous) models to be attractive like even with 
large z values, hence it needs to be investigated. The situation, however, is not necessarily as bad as it looks, since 
even though pressure at first seems to make the situation worse, it may actually recover it. This is due to the chances 
pressure brings with it to the evolution equation. For example, if the evolution equation changes from second order 
differential equation to higher order one, it may chance the structure of the stability of the system dramatically. This 
may even happen for homogeneous models, thus the result here received for dust and dark energy filled universe 
should not be interpreted as final. 
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